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Abstract 

In this paper we show how to simulate and estimate a COGARCH(p,q) model in the 
R package yuima. Several routines for simulation and estimation are available. Indeed for 
the generation of a COGARCH(p,q) trajectory, the user can choose between two alternative 
schemes. The first is based on the Euler discretization of the stochastic differential equations 
that identifies a COGARCH(p,q) model while the second one considers the explicit solution 
of the variance process. 

Estimation is based on the matching of the empirical with the theoretical autocorrelation 
function. In this case three different approaches are implemented: minimization of the mean 
square error, minimization of the absolute mean error and the generalized method of moments 
where the weighting matrix is continuously updated. 

Numerical examples are given in order to explain methods and classes used in the yuima 
package. 
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1 Introduction 

The Continuous-Time GARCH(1,1) process has been introduced in |18| as a continuous coun¬ 
terpart of the discrete-time GARCH (1,1) model proposed by [3]. 

The idea is to develop in continuous time a model that is able to capture some stylized facts 
observed in financial time series m exploiting only one source of randomness for returns and 
for variance dynamics. Indeed, in the Continuous-Time GARCH (COGARCH hereafter), the 
stochastic differential equation for variance is driven by the discrete part of the quadratic 
variation of the same Levy process used for modeling returns. The continuous nature of the 
COGARCH makes it particularly appealing for discribing the behaviour of high frequency 
data [see m for an application of method of moments using intraday returns]. 

The generalization to higher order COGARCH(p,q) processes has been proposed in 0[S]. 
Starting from the observation that the variance of a GARCH(p,q) is an ARMA(q, p-1), the 
Variance is modeled with a CARMA(q,p-l) process [see [7l 1271 16| and many others] driven by 
the discrete part of the quadratic variation of the Levy process in the returns. Although this 
representation is different from the one used by [18] for the COGARCH(l,l) process, this last 
can be again retrieved as a special case. 

Many authors recently have investigated the COGARCH(l,l) model from a theoretical and 
an empirical point of view [see |22[ 116112^ 1^ and many others]. Some R codes for estimation 
and simulation of a COGARCH(l,l) driven by a Compound Poisson and Variance Gamma are 
available in m- For the general COGARCH(p,q), the main contribution remain the seminal 
works [5] and [9]. The aim of this paper is to describe the simulation and the estimation 
schemes in the yuima package |26| for a GOGARGH(p,q) model driven by a general Levy pro¬ 
cess. Based on our knowledge yuima is the first R package available on CRAN that allows the 
user to manage a higher order GOGARGH(p,q) model. Moreover, the estimation algorithm 
gives the option to recover the increments of the underlying noise process and estimates the 
Levy measure parameters. We recall that a similar procedure is available in yuima also for 
the GARMA model ]?, see]for a complete discussion]IacusMercur2015. The yuima package is 
developed within the YUIMA project [8] whose aim is to provide an environment for simulation 
and estimation of stochastic differential equations. 

The outline of the paper is as following. In Sect. [2] we discuss the main properties of the 
COGARCH(p,q) process. In particular we review the condition for existence of a strictly 
stationary variance process, its higher moments and the behaviour of the autocorrelation of 
the square increments of the COGARCH(p,q) model. In Sect. [3] we analyze two different 
simulation schemes. The first is based on the Euler discretization while the second one uses 
the solution of the state process in the CARMA(q,p-l) model. Sect. [I] is devoted to the esti¬ 
mation algorithm. In Sect. [5] we show the main classes and corresponding methods in yuima 
package and in the Sect, [^we present some numerical examples about the simulation and the 
estimation of a COGARCH(p,q) model 
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2 COGARCH Models driven by a Levy process 

In this section we review the mathematical definition of a COGARCH(p,q) process and its 
properties. In particular we focus on the conditions for the existence of a strictly stationary 
COGARCH(p,q) process and compute the first four unconditional moments. The existence of 
higher order moments plays a central role for the computation of the autocorrelation function 
of the squared increments of the COGARCH(p,q) model and consequently the estimation 
procedure implemented in the yuima package. 

The COGARCH(p,q) process, introduced in [Q as a generalization of the COGARCH(l,l) 
model, is defined through the following system of stochastic differential equations: 

r dGt^VVtdLt 

I Vt = ao + (1) 

[ dYt = AYt-dt + e{ao + aJYt-)d[L,L]f 

where q and p are integers such that q > p > 1. The state space process Yt is a vector with q 
components: 

Yt = [yi,t,...,y„,t]^. 

The vector a £ TV^ is defined as: 


— [ai,..., CLp^ Up+i,..., 


with Op+i = ■ • ■ = Oq = 0. The companion q x q matrix A is 

■ 0 1 ... 0 

A= : : : 

0 0 ... 1 

-bg —bq-1 ... —6i 

The vector e £ TV^ contains zero entries except the last component that is equal to one. 

is the discrete part of the quadratic variation of the underlying Levy process Lt and is 
defined as: 

[L,L]f:= ^ (AL.)^ (2) 

0<s<t 


Remark 1 A COGARCH(p,q) model is constructed starting from the observation that in the 
GARCH(p,q) process, its discrete counterpart, the dynamics of the variance is a predictable 
ARMA(q, p-1) process driven by the squares of the past increments. In the COGARCH(p,q) 
case, the ARM A process leaves the place to a CARMA(q,p-l) model [see for details about 
the CARMA(p,q) driven by a Levy process] and the infinitesimal increments of the COCA- 
RCH(p,q) are defined through the differential of the driven Levy Lt as done in ([T|l. 

As observed above the COGARCH(p,q) model generalizes the COGARCH(l,l) process that 
has been introduced following different arguments from those for the (p, q) case. However 
choosing q = 1 and p = 1 in (O the GOGARGH(l,l) process developed in (TH] US] can be 
retrieved through straightforward manipulations and, for obtaining the same parametrization 
in Proposition 3.2 of [18], the following equalities are necessary: 

uJo = aobi, LOi = and 77 = 61. 

Before introducing the conditions for strict stationarity and the existence of unconditional 
higher moments, it is worth noting that the state space process Yt can be seen as a Multi¬ 
variate Stochastic Recurrence Equation and the corresponding theory can be applied to the 
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COGARCH(p,q) process [see 0 more details] in order to derive its main feature^. In the 
case of the Compound Poisson driven noise, the representation through the stochastic differ¬ 
ence equations is direct in the sense that the random coefficients of the state process Yt can 
be written explicitly while in the general case, it is always possible to identify a sequence of 
Compound Poisson processes that as limit to the choosen driven Levy process. 

In the following, we require that matrix A can be diagonalized, that means: 

A = SDS~^ 


with 



■ 1 

1 ■ 


Ai 

s = 

Ai 

. Aq 



. Ar^ . 

. Ar'. 


Aq . 


where the Ai, Aq,..., A, are the eigenvalues of matrix A and are ordered as follows: 


(3) 


K{Ai} > K{A2} > ... > K{Aq}. 

Applying the theory of stochastic recurrence equations, provide A sufficient condition for 
the strict stationarity of a COCARCH(p,q) model. We review the result for the stationarity of 
the state process Yt. Two fundamental assumptions are the fact that the eigenvalues Ai,..., Aq 
are distinct and the underlying process L must have a non-trivial i/l (1) measure. Then, the 
process Yt converges in distribution to the random variable Yoo if exist some r £ [1, -l-oo] such 
that: 

/ + 00 

ln(l-b ||S-^ea^S'||,l^)d!zn ( 1 ) < K{Ai} (4) 

-OO 

for some matrix S such that the matrix A is diagonalizable. If we choose as a starting condition 
Yo = Yoo than the process Yt is strictly stationary and consequently the variance Vt is also 
a strictly stationary process. Unlikely for the general case, the inequality in gives only a 
sufficient condition about the strict stationarity, and it is difficult to verify in practic^EI. As 
shown in m and remarked in , the condition in @ is also necessary for the COCARCH(l,l) 
case and can be simplified as: 



-I- mP) duL (I) < bi. 


(5) 


Using again the SRE theory, it is also possible to determine the condition for existence of higher 
moments of the state process Yt. In this way it is possible to determine the autocorrelations 
of the squared COCARCH increments that are used in the estimation procedure illustrated 
in Section 21 As reported in [^, the tc-th order moment of the process Yt exists finite if, for 
some r > 1 and S such that the matrix A is diagonalizable, the following conditions hold: 

<+oo, r [(l-b||S-^ea^S||,l^)''-l]diZi(/)<K{Ai}tv. (6) 


^The Stochastic Recurrence Equations theory HE] has been also used to prove the strictly and weakly station¬ 
arity for the GARCH(p,q) model ^ 

^In the yuima package a diagnostic for condition 0 is available choosing matrix S as done in 0 and r = 2. We 
remark that the process Yt is stationary if the diagnostic gives a positive answer otherwise we can conclude nothing 
about the stationarity of the process. 
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As special case of the unconditional stationary mean E (yoo) = —aofj, {A + fieaJ) ^ e of 
the vector process It exists if 

E {lI) < +0O, ||S'-^ea^5||rAt < K{Ai} 


where 


M := 



l^dVL (1) 


(7) 


is the second moment of the Levy measure vl (I). It is worth noting that the condition ^ 
ensures also the strict stationarity since, using ln(l + a;) < a;, we have the following sequences 
of inequalities: 



ln(l + ||S ^ea’'’S||r^^) dia^ (/) < 115 ^ea'''5||r/a < K{Ai} . 


For the stationary covariance matrix [9] 
cov (Voo) = 


O-obqP 


■f 


At T 

e ee'e at, 


{bq — fiaif (1 — m) 
the existence of the second moment of Yoo in m becomes: 

E (Lf) < oo, ||5"^ea^5||,p < 2 (-K{Ai} - ||5“^ea^5||,M) 


( 8 ) 


where 

/ + 00 

l^di^L (0 

- OO 

is the fourth moment of the Levy measure (1) and 

T At T t j, 

a'e ee'e adt. 

Before introducing the higher moments and the autocorrelations, we recall the conditions for 
the nonegativity for a strictly stationary variance process in ©■ Indeed, under the assumption 
that all the eigenvalues of matrix A are distinct and the relation in @ holds, the variance 
process Vt > ao > 0 a.s. if: 

>0, Vt > 0. (9) 

The condition in is costly since we need to check it each time. Nevertheless some useful 
controls are available |28| . 

1. A necessary and sufficient condition to guarantee that > 0 in the COGARCH(2,2) 

case is that the eigenvalues of A are real and 02 > 0 and a\ > —a 2 \ (A) where A (A) is 
the biggest eigenvalue. 

2. Under condition 2 < p < q, that all eigenvalues of A are negative and ordered in an 
increasing way Ai > A 2 Ap_i and 7 j are the roots of a ( 2 ) = 0 ordered as 

0 > 71 > 72 > ... > 7 p-i. Then sufficient condition for ® is 

k k 

Vi < VfcG {!,...,p-1}. 

3. For a COGARCH(l,q) model a sufficient condition that ensures (|9l) is that all eigenvalues 
must be real and negative. 
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Combining the requirement in Q with that in @ it is possible to derive the higher moments 
and the autocorrelations for a COGARCH(p,q) model. As a first step, we dehne the returns 
of a COGARCH(p,q) process on the interval + r], Vt > 0 as: 

( 10 ) 

Let Lt be a symmetric and centered Levy process such that the fourth moment of the associated 
Levy measure is finite, we define the matrix A as: 

A := A + fiea^ . 


It is worth noting that the structure of A is the same of A except for the last row q where 
Aq^ = {—bq + /iui,..., —b\ + /ittq). For any t > 0 and for any h > r > 0, the first two moments 
of (fTUl) are 



( 11 ) 


and 


E 


(Gr^y]=j^^E[Ll]. ( 12 ) 

L V / J bq — fiai '■ ■' 

The computation of the autocovariances and variance for squared returns Gnu require the 
following quantities defined as: 


Po = [si-i (A-^ - r/) - /] cov (coo) 

Ph = (/ - A-i - /j cov (e^) ’ 

Qo = 6m [(r7 - - l)') cov (coc) - (A’I - /) - r/) cov (c^) i"] e 

Qh = (/ - [(/ - - A-^ - ij cov (e^) A"] e 

(14) 

and 

R = 2r^ pA + p. (15) 

The terms Ph and Po are g x g matrices. Qh and Qo are g x 1 vectors and the term P is a 
scalar. The g x g matrix cov (coo) is defined as: 


cov 



At T 1, 

e ee' e at 


(16) 


where p, p and m have been defined before. 

Using the g x g matrix Ph and the g x 1 vector Qh the autocovariances of the squared returns 
uni are defined as: 


flofeg (a'^P^a + aQh) , ^ 
-^-2-. 


7 r (h) ~ cov [(Gt’’^)^, (G(V)^1 = 
^ ^ t t , V t+hl J 

while the variance of the process is 

7.(0):=.a.((G, )j= 


(17) 


(18) 
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Combining the autocovariances in (nzi) with m and m we obtain the autocorrelations: 


(h) ■= > (^) _ (a^Pfea + aQft) 

’ Jr (0) (aT Poa + aQo + R) 


(19) 


We conclude this Section with the following remark. As done for the GARCH(p,q) model the 

fr") 

autocovariances of the returns G) are all zeros. 


3 Simulation of a COGARCH(P,Q) model 

In this Section we illnstrate the theory behind the simulation routines available in the yuima 
package when a COGARCH(p,q) model is considered. The corresponding sample paths are 
generated according two different schemes. 

The first method, in the Yuima project [S], is based on the Euler-Maruyama m discretization 
of the system in ©• In this case the algorithm follows these steps: 

1. We determine an equally spaced grid from 0 to T composed by N intervals with the same 
length At 

0,..., (n — 1) X At, n x At,... ,N x At. 

2. We generate the trajectory of the nnderlying Levy process Lt sampled on the grid defined 
in step 1. 

3. Choosing a starting point for the state process Yb = yo and Go = 0, we have 

= (/ + AAt) y„_i + e (ao + a^y„_i) A [LL]^ . (20) 

The discrete quadratic variation of the Lt process is approximated by 

A[LL]^ = (L„-L„_l)^ (21) 

4. Once the approximated state process Yn is obtained we can generate the trajectory of 
Variance and the process Gn according the following equations: 

Vn = O>0 “t“ ^"^Yn—l (22) 

and 

G„ = G„-1 + (in - Tn-i). (23) 

It is worth noting that, although the discretized version of the state process Yn in (Uni) can 
be seen as a stochastic recurrence equation, the conditions for stationarity and non-negativity 
for the variance V„ process are not the same of ones analyzed in the Section [21 In particular it 
is possible to determine an example where the discretized variance process Vn assumes negative 
values while the true process is non negative with probability one. 

In order to clarify deeper this issue we consider a COGARCH(l,l) model driven by a Variance 
Gamma Levy process [see |21[ I19| for more details about the VG model]. In this case, the 
condition for the non-negativity of the Variance in m is ensured if uq > 0 and ai > 0 while 
the strict stationarity condition in (O for the COGARCH(l,l) is E[L‘^] — 1 and a\ — b\ < 0. 
The last two requirements guarantee also the existence of the stationary unconditional mean 
for the process Vt. We define the model using the yuima function setCogarch. Its usage is 
completely explained in the Section [S] The following command line instructs yuima to build 
the COGARCH(l,l) model with VG noise: 
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> modell <- setCogarchCp = 1, q = 1, work = FALSE, 

+ measure=list("rngaimiia(z, 1, sqrt(2), 0, 0)"), measure.type = "code", 

+ Cogarch.var = "G", V.var = "v", Latent.var="x", XinExpr=TRUE) 

Choosing the following values for the parameters 

> paraml <- listCal = 0.038, bl = 301, aO =0.01, xOl = 0) 

the COGARCH(l,l) is stationary and the variance is strictly positive. Nevertheless, if we 
simulate the trajectory using the Euler discretization, the value of At can lead to negative 
values for the process as shown in the example; 

> Terminall=5 

> nl=750 

> sampl <- setSampling(Terminal=Terminall, n=nl) 

> set.seed(123) 

> siml <- simulate(modell, sampling = sampl, true.parameter = paraml, 

+ method="euler") 

> plot(siml, main="Sample Path of a VG C0GARCH(1,1) model with Euler scheme") 


Sample Path of a VG COGARCH(1,1) model with Euler scheme 



t 

Looking to the Figure, we observe a divergent and oscillatory behaviour for the simulated 
state Yn and the variance V„ processes while, from theoretical point of view, the conditions 
for nonnegativity and stationarity for the variance of a COGARCH(l,l) model are satisfied 
by the values used for the parameters. 

To overcome this problem, we provide an alternative simulation scheme based on the solution 
of the process Yn given the starting point Yn-i- 

Applying the Ito’s Lemma for semimartingales [25] to the transformation e~^*Yt, we have: 

= Ft-At-/ Ae-^"Yn-du+ [ e"^“dy„+^ (F, - F.-) - (F, - F,- 

Jt-At Jt-At ^ 

We substitute the definition of Yt in m and get 

e-^*F = F_At-/ Ae-^'^Yu-Au+[ e"^“AF„_du+/" (uo + a'F„-) d 

Jt-At Jt-At Jt-At 
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Using the following property for an exponential matrix 

= A ^/ ++ ...^ 

= (^A + tA^+ 

= ^7 + t^++ ...^ ^ = 

we get 

Yt = e^^Yt-At+ f e'^(*-“^e(ao + a'Uu_)d[LL];^ (24) 

Jt-At 

Except for the case where the noise is a Compound Poisson, the simulation scheme follows 
the same steps of the Euler-Maruyama discretization where the state space process Yn on the 
sample grid is generated according the approximation of the relation in (|24l): 

K = (no + n'K-i) (25) 

or equivalently: 

y„ = [LL]^ + (/ + ea'A [LL]^) y„_i. (26) 

where A := [LL]'^ — [LL]'^_^ is the increment of the discrete part of quadratic variation. 

In the previous example, the sample path is simulated according the recursion in (EH) choosing 
method = mixed in the simulate function as done below: 

> set.seed(123) 

> sim2 <- simulate(modell, sampling = sampl, true.parameter = paraml, 

+ method="mixed") 

> plot(sim2, main="Sample Path of a VG CDGARCH(1,1) model with mixed scheme") 

In the case of the COGARCH(p,q) driven by a Compound Poisson Levy process, a trajec¬ 
tory can be simulated without any approximation of the solution in (EH). Indeed it is possible 
to determine the time where the jumps occu]0 and then evaluate the corresponding quadratic 
variation in an exact way. Once the trajectory of a random time is obtained, the piecewise 
constant interpolation is used on the fixed grid, in order to mantain the cadlag property of 
the sample path. 


4 Estimation of a COGARCH(P,Q) model in the 
yuima package 

In this Section we explain the estimation procedure that we propose in the yuima package for 
the COGARCH(p,q) model. As done for the CARMA(p,q) model driven by Levy process |15 |. 
even in this case a three step estimation procedure is proposed that allows the user to obtain 
estimated values for the GOGARGH(p,q) and the parameters of the Levy measure. 

This procedure is structured as follows: 

®In a general Compound Poisson the jump time follow an exponential r.v. with rate A 
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Sample Path of a VG C0GARCH(1,1) model with mixed scheme 



t 


1. Using the moment matching procedure explained below, we estimate the COGARCH(p,q) 
parameters a [ai,..., a^], b ;= [6i, ... ,bq] and the constant term ao in the Variance 
process Vt. In this phase, the estimation is obtained by minimizing some distances be¬ 
tween the empirical and theoretical autocorrelation function. 

2. Once the COGARCH(p,q) parameters are available, we recover the increments of the 
underlying Levy process using the methodology describe below. 

3. In the third step, we use the increments obtained in the second step and estimate the 
Levy measure parameters by means of maximum likelihood estimation procedure. 


Let Go, Gi,..., G„, ..., Gt be the observed values of the process subsampled at equally 
spaced instants 0, h, 2h ,..., Nh where the length of each interval is /i = The time of the 
last observation is T and N is the number of observations of the process Gn- 
In the following we assume that the underlying Levy process is symmetric and centered in 
zero. Starting from the sample {Gnj^^Q we define the GOGARGH(p,q) increments of lag one 
as: 


Gi^^ ~ G„ - Gr, 


and the increments of lag r as: 

G^ := G„ - Gr^-r (27) 

where r > 1 is a natural number and for r = 1 the definition in (H) coincides with 113. 

(r) 

It is worth mentioning that the increments Gn can be obtained as a time aggregation of 
increments of lag one as follows: 

n 

gM= G« (28) 

h=n — r 

The time aggregation in ([281) can be useful during the optimization routine when the values 
of increments G^^^ are very small in absolute value. 
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we compute the empirical second moment 


Using the 


sample 


N — d — r + 


iv —a 


and empirical autocovariance function 7 {h) is defined as: 

^ N — d / 2 \ 

7 r- {h) ■.= ^ ^ r \ 1 ^ V “Ar) h^0,l,...,d 

n=r ' ^ 

where d is the maximum lag considered. 

The empirical autocorrelations are: 


pT Qi) 


It jh) 

7i" ( 0 ) 


(29) 


We use the theoretical and empirical autocorrelations in order to define the moment conditions 
for the estimation procedure. By the introduction of the q + p vector 0 := (a, b) we define the 
vector function g {6) : R'^ as follows: 


g{e) :=£;[/ 


(30) 


where / (^Gn'^ , a, b^ is a d dimensional real function where the components are defined as: 


fH[AGl^\e^=Pr (h) 


-PtJUgG^) 


jll 


It ( 0 ) 


h = l,...,d. (31) 


In the estimation algorithm, we replace the expectation in with the sample counterpart. 
The components of vector g {9) = [gi (d) ,..., ga (d)]’^ are defined as: 


9h (9) = 


N-d 

l-r+l ^ 


N-d-' 

= Pr{h) - Pr{h), h= 1,..., d. 


Pr (h) - 


{Gi^y-Pr^ ((gW)'-A 

7r (0) 


(32) 


The vector 9 — (a, b) containing the COGARCH(p,q) parameters are obtained by minimizing 
some distances between empirical and theoretical autocorrelations. The optimization problem 
is: 

min d{pr,Pr) 
e^mi+p 

where the components of vectors pr and pr are the theoretical and empirical autocorrelations 
defined in m and 1291) respectively. Function d {x, y) measures the distance between vectors 
X and y. In the yuima environment, three distances are available and listed below: 

1 . the Li norm 

IIA(0)l|i=El5'*WI- (33) 

h = l 
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2. the squared of L 2 norm 


d 

= (34) 

h=l 

3. The quadratic form 

ll5Wllw=5W^W5(0) (35) 

where the positive definite weighting matrix W is choosen to obtain efficient estimators 
between those that belong to the class of asymptotically normal estimators. 

It is worth noting that the objective function \\g (6')||| is a special case of the function \\g (0)||w 
where the weighting matrix W coincides with the identity matrix. Both distances are related 
with the Generalized Method of Moments (GMM) introduced by |11| . Under some regularity 
conditions |24| . the GMM estimators are consistent and, for any general positive definite 
matrix W, their asymptotic variance-covariance matrix V are: 

V = —--(D^WD)"^ D^WSWD (D^WD)”^, 

N-d-r +1 ^ ^ ’ 


The matrix D is defined as: 


D = E 


df (G(r\0) 


While 

S = E[f (gW,0)/(gM. 0)'] . 
For the squared L 2 norm in (1341) matrix V becomes: 


(36) 

(37) 


V = 


N — d — r + 1 


(D^D)~^D^SD(D^D) 


(38) 


while for (13511 . as observed above, the choice of the matrix W is done in order to obtain 
efficient estimators in the class of all asymptotically normal estimators. To obtain this result 
we prefer to use the Continuously Updated GMM estimator m In this case the matrix W is 
determinated simultaneusly with the estimation of the vector parameters 9. Introducing the 
function ||g (^)||^ as the sample counterpart of the quadratic form in (1351) . the minimization 
problem becomes: 

min ||<7(e)||^=fl(0)^W(e)Ke) 

06R9+P 

where function W (9) maps from 77.^+'^ to and is defined as: 

^ N—r—d 


W(6l) = 


N — r — d + 1 


E f {Gi:\e) f {Gi^\9y 


(39) 


Observe that W (9) is a consistent estimator of matrix S ^ that means: 


W{9) S"" 

N-t+oo 

consequently the asymptotic variance-covariance matrix V in (1381) becomes: 


(40) 


V = 


N — r — d + 1 


(D^S'^D) 


(41) 


Once the estimates of vector 9 are obtained, the user is allowed to retrieve the increments 
of the underlying Levy process according the following procedure. This stage is independent 
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on the nature of the Levy measure but it is only based on the solution of the state process Yt 
and on the approximation of the quadratic variation with the squared increments of the Levy 
driven process. 

Starting from the discrete time equally spaced observations Gi, G 2 ,. •., GvAt, we remark 
that the increment AG* := Gt — Gt- can be approximated using the observations {G„At}^^o 
as follows: 

AGt ~ AGnAt = GnAt — G(„_l)At (42) 

Recalling that Gt = X)o<s<t i^Ls), the approximation in (BH) becomes: 

AGt « \/V^t {AL nAt ) (43) 

where V„At is the value of the variance process at the observation time t = and ALnAt = 
L„At — L(n-i)At is the discretized Levy increment from t — At to t. 

Using the discretization scheme introduced in (1251) the process Y^st is written as: 

Yr^st ~ (ao + aUt-At) ([LL]" - [LL](„_i)) (44) 

since the difference [LL]"^ — ~ (AL„At)^ and using the result in (14311 . the difference 

equation dm can be approximated in terms of the squared increments of COGARCH(p,q) 
process and we have: 

Yt « e^^^Yt-At + (ao + a'Yt-At) (ALtf .... 

= e^^V,_At+e^(^*)e(AGO^ 


Choosing Vo equal to the unconditional mean of the process Yt, we are able to retrieve 
its sample path according to the recursive equation in (H51) . The only quantities that we 
need to compute are the squared increments of the COGARCH(p,q) process on the grid 
{0, At, 2At,..., nAt..., NAt}. The estimated state process in (1451) is also useful for getting 
the estimated trajectory of the variance process. Finally note the Levy increment at a general 
time t = nAt is obtained as: 


AIjn5t — 


AGnSt 

^Yn6t 


(46) 


Once the increments of the underlying Levy are obtained, it is possible to obtain the estimates 
for the Levy measure parameters through the Maximum Likelihood Estimation procedure [we 
refere to the yuima documentation [8] [26] for the available Levy processes]. 


5 Package R 

In this Section we illustrate the main classes and methods in yuima that allow the user to deal 
with a COGARCH(p,q) model. The routines implemented are based on the results considered 
in Section [3] for simulation and Section [4] for the estimation procedures. 

In particular we classify these classes and methods in three groups. The first group contains 
the classes and functions that allow the user to define a specific COGARCH(p,q) model in 
the yuima framework. The second group is used for the simulation of the sample paths for the 
COGARCH(p,q) model and the third is related to the possibility of estimation using simulated 
or real data. 
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5.1 Classes and Methods for the definition of a COGARCH(P,Q) 
model 

The main object for a COGARCH(p,q) process in the yuima environment is an object of 
yuima. cogarch-class that contains all the information about a specific COGARCH(p,q) pro¬ 
cess. This class extends the yuima.model-class and it has only one additional slot, called 
Oinfo, that contains an object of cogarch. info-class. We refer to the yuima documentation 
for a complete description of the slots that constitnte the objects of class yuima.model. In 
this paper we focus only on the object of class cogarch.info. In particular its structure is 
composed by the slots listed below: 

• @p is an integer that is the number of moving average coefficients in the Variance process 

II. 

• @q is an integer number that corresponds to the dimension of autoregressive coefficients 
in the variance process Vt- 

• @ar. par contains a string that is the label for the autoregressive coefficients. 

• Sma.par is the Label for the moving average coefficients. 

• Sloc.par indicates the name of the location coefficient in the process Vt. 

• SCogarch.var string that contains the name of the observed process Gt- 

• @V . var is the Label of the Vt process. 

• SLatent.var indicates the label of the state process Yt- 

• SXinExpr is a logical variable. If the value is FALSE, default value, the starting condition 
for the state process Yt is a zero vector. Otherwise the user has to fix the starting 
condition as argument true.parameter in the method simulate. 

• ^measure identifies the Levy measure of the underlying noise and consequently the dis¬ 
crete part of the quadratic variation that drives the state process. 

• ^measure. type says if the Levy measure belongs to the family of Compound Poisson or 
is another type of Levy 

The user builds an object of class yuima. cogarch through to the constructor setCogarch: 

setCogarchCp, q, ar.par = "b", ma.par = "a", loc.par = "aO", Cogarch.var = "g", 
V.var = "v", Latent.var = "y", jump.variable = "z", time.variable = "t", measure 
= NULL, measure.type = NULL, XinExpr = FALSE, startCogarch = 0, work = FALSE, ...) 

The arguments used in a call of the function setCogarch are illustrated in the following 

list: 

• p: A non negative integer that is the number of the moving average coefficients in the 
variance process. 

• q: A non negative integer that indicates the number of the autoregressive coefficients in 
the variance process. 

• ar.par: A character-string that is the label of the autoregressive coefficients. 

• ma.par: A character-string that is the label of the autoregressive coefficients. 

• loc.par: A string that indicates the label of the location coefficient in the variance 
process. 

• Cogarch. var: A character-string that is the label of the observed COGARCH process. 
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• V.var: A character-string that is the label of the latent variance process. 

• Latent.var: A character-string that is the label of the latent process in the state space 
representation for the variance process. 

• jump. variable: Label of the underlying Levy process . 

• time. variable: Label of the time variable. 

• measure: Levy measure of the underlying Levy process. 

• measure. type: Label that indicates whether the underlying noise is a Compound Poisson 
process or another Levy without the diffusion component. 

• XinExpr: A logical variable that identifies the starting condition. In particular, the 
default value XinExpr = FALSE implies that the starting condition for the state process 
is zero. Alternatively XinExpr = TRUE means that the user is allowed to specify as 
parameters the starting values for each component of the state variable. 

• startCogarch: Initial value for the COGARCH process. 

• ...: Arguments to be passed to setCogarch such as the slots of the yuima.model-class. 

5.2 Classes and Methods for the simulation of a COGARCH(P,Q) 
model 

The simulate is a method for an object of class yuima.model. It is also available for an object 
of class yuima. cogarch. The function requires the following inputs: 

simulate(object, nsim=l, seed=NULL, xinit, true.parameter, space.discretized 
= FALSE, increment.W = NULL, increment.L = NULL, method = "euler", hurst, methodfGn 
= "WoodChan", sampling=sampling, subsampling=subsampling, ...) 

In this work we focus on the argument method that identifies the type of discretization 
scheme for the time when the object belongs to the class of yuima.cogarch. The default 
value euler means that the simulation of a sample path is based on the Euler-Maruyama 
discretization of the stochastic differential equations. This approach is available for all objects 
of class yuima.model. For the COGARCH(p,q) an alternative simulation scheme is available 
choosing method = mixed. In this case the generation of trajectories is based on the solution 
(IMl) for the state process. In particular if the underlying noise is a Gompound Poisson Levy 
process, the trajectory is built using a two step algorithm. First the jump time is simulated 
internally using the Exponential distribution with parameter A and then the size of jump is 
simulated using the random variable specified in the slot yuima. cogarchOmodelSmeasure. For 
the other Levy processes, the simulation scheme is based on the discretization of the state 
process solution (l25ll in Section [5] 

5.3 Classes and Methods for the estimation of a COGARCH(P,Q) 
model 

The cogarch.gmm class is a class of the yuima package that contains estimated parameters 
obtained by the gmm function. 

• Omodel is an object of of yuima. cogarch-class. 

• OobjFun is an object of class character that indicates the objective function used in 
the minimization problem. L2 refers to the squared of L 2 norm in (1 34 II . L2CUE for the 
quadratic form (I35II and the LI for the Li norm in (I33II 

• @call is an object of class language. 
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• Scoef is an object of class numeric that contains the estimated parameters. 

• Sfullcoef is an object of class numeric that contains the estimated and hxed parameters 
from the user. 

• @vcov is an object of class matrix. 

• @min is an object of class numeric. 

• Sminuslogl is an object of class function. 

• ^method is an object of class character. 

The cogarch. gmm. incr is a class of the yuima package that extends the cogarch. gmm-class 
and is Hlled from the function gmm. 

• Incr.Lev is an object of class zoo that contains the estimated increments of the noise 
obtained using cogarchNoise. 

• modelis an object of yuima. cogarch-class. 

• logL.Incr is an object of class numeric that contains the value of the log-likelihood for 
the estimated Levy increments. 

• objFun is an object of class character that indicates the objective function used in the 
minimization problem. The values are the same for the slot OobjFun in an object of class 
cogarch.gmm. 

• call is an object of class language. 

• coef is an object of class numeric that contains the estimated parameters. 

• fullcoef is an object of class numeric that contains estimated and fixed parameters. 

• vcov is an object of class matrix. 

• min is an object of class numeric. 

• minuslogl is an object of class function. 

• method is an object of class character. 

The function gmm returns the estimated parameters of a COGARCH(p,q) model. The pa¬ 
rameters are obtained by matching the theoretical with the empirical autocorrelation function. 

gmmCyuima, data = NULL, start, method="BFGS", fixed = listO, lower, upper, lag.max 
= NULL, equally.spaced = TRUE, Est.Incr = "Noincr", objFun = "L2") 

• yuima is a yuima object or an object of yuima. cogarch-class 

• data is an object of class yuima.data-class contains the observations available at uni¬ 
formly spaced instants of time. If data=NULL, the default, the function uses the data in 
an object of yuima-class. 

• start is a list that contains the starting values for the optimization routine. 

• method is a string that indicates one of the methods available in the function optim. 

• fixed a list of fixed parameters in the optimization routine. 

• lower a named list for specifying lower bounds for parameters. 

• upper a named list for specifying upper bounds for parameters. 

• lag.max maximum lag for which we calculate the theoretical and empirical acf. Default 
is y/N where N is the number of observations. 


16 



• equally. spaced Logical variable. If equally. spaced = TRUE, in each unitary interval 
we have the some number of observations. If equally. spaced = FALSE, each unitary 
interval is composed by different number of observations. 

• Est. Incr a string variable. If Est. Incr = "Noincr ", default value, girnn returns an object 
of class cogarch.gimn-class that contains the COGARCH parameters. If Est.Incr = 
"Incr" or Est. Incr = "IncrPar" the output is an object of class cogarch. gimn. incr-class. 
In the first case the object contains the increments of the underlying noise while in the 
second case it contains also the estimated parameters of the Levy measure. 

• objFun a string variable that indentifies the objective function in the optimization step. 
objFun = "L2", default value, the objective function is a quadratic form where the 
weighting matrix is the identity one. objFun = "L2CUE" the weighting matrix is es¬ 
timated using Continuously Updating GMM (L2CUE). objFun = "LI", the objective 
function is the mean absolute error. In the last case standard errors for estimators are 
not available. 

Function gmm uses function cogarchNoise for the estimation of the underlying Levy in a 
COGARCH(p,q) model. This function assumes that the underlying Levy process is symmetric 
and centered in zero. 

cogarchNoise(yuima.cogarch, data=NULL, param, mu=l) 

The arguments of the cogarchNoise are listed below 

• yuima.cogarch is an object of yuima-class or an object of yuima.cogarch-class that 
contains all the information about the COGARCH(p,q) model. 

• data is an object of class yuima.data-class that contains the observations available at 
uniformly spaced instants of time. If data=NULL, the default, the cogarchNoise uses the 
data in an object of yuima.data-class. 

• param is a list of parameters for the GOGARCH(p,q) model. 

• mu is a numeric object that contains the value of the second moments of the Levy measure. 

6 Numerical results 

In this section we show how to use the yuima package for the simulation and the estimation of a 
GOGARGH(p,q) model driven by different symmetric Levy processes. As a first step we focus 
on a GOGARGH(l,l) model driven by different Levy processes available on the package. In 
particular we consider the cases in which the driven noise are a Gompound Poisson with jump 
size normally distributed and a Variance Gamma process. In the last part of this section, we 
show also that the estimation procedure implemented seems to be adequate even for higher 
order GOGARGH models. In particular we simulate and then estimate different kind of 
GOGARGH(p,q) models driven by a Compound Poisson process where the distribution of the 
jump size is a normal. 

6.1 Simulation and Estimation of a COGARCH(l,l) 

The first example is a COGARCH(l,l) model driven by a Compound Poisson process. As a 
first step, we choose the set of the model parameters: 

> numSeed <- 200 

> param.cp <- listCal = 0.038, bl = 0.053, aO = 0.04/0.053, 

+ lambda = 1, eta=0, sig2=l, xOl = 50.33) 
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ai, bi and oq are the parameters of the state process Yt. A is the intensity of the Compound 
Poisson process while rj and are the mean and the variance of the jump size. xo,i is the 
starting point of the process Xt, the choosen value is the stationary mean of the state process 
and it is used in the simulation algorithm. 

In the following command line we define the model using the setCogarch function. 

> mod.cp <- setCogarch(p = 1, q = 1, work = FALSE, 

+ measure=list (intensity="laiiibda" ,df=list ("dnormCz ,eta, sig2) ")) .measure .type = "CP", 

+ Cogarch.var = "g", V.var = "v", Latent.var="x", 

+ XinExpr=TRUE) 

We simulate a sample path of the model using the Euler discretization. We fix At = and 
the command lines below are used to instruct yuima for the choice of the simulation scheme: 

> Term <- 1600 

> num <- 24000 

> set.seed(numSeed) 

> samp.cp <- setSampling(Terminal=Term, n=num) 

> sim.cp <- simulate(mod.cp, true.parameter=param.cp, 

+ sampling=samp.cp, method="euler") 

In the following figure we show the behaviour of the simulated trajectories for the COGA- 
RCH(1,1) model Gt, the variance Vt and the state space Xt'. 

> plot(sim.cp, main = "simulated C0GARCH(1,1) model driven by a Compound Poisson process") 



We use the two step algorithm developed in Section |4] for the estimation of the COG- 
ARCH(p,q) and the Levy measure parameters. In the yuima function gmm, we fix objFun 
= L2 meaning that the objective function used in the minimization is the mean squared er¬ 
ror. Setting also Est. Incr=IncrPar, the function gimn returns the estimated increments of the 
underlying noise. 

> res.cp <- gimn(sim.cp, start = param.cp, objFun = "L2", Est. Incr = "IncrPar") 

The results can be displayed using the method summary and in the following figure we report 
the recovered increments of the underlying noice process. 

> summary(res.cp) 

Two Stages GMM estimation 
Call: 

gmm(yuima = sim.cp, start = param.cp, Est.Incr = "IncrPar", objFun = "L2") 
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Coefficients: 


bl 

Estimate 

6.783324e-02 

Std. Error 

0.06862392 

al 

3.403071e-02 

0.02897625 

aO 

1.032014e+00 

NA 

lambda 

1.073912e+00 

NA 

eta 

6.818470e-09 

NA 

sig2 

7.837838e-01 

NA 


Log.objFun L2: -3.491179 


Number of increments: 24000 
Average of increments: -0.002114 
Standard Dev. of increments: 0.256610 


-2 log L of increments: 2851.529874 

Summary statistics for increments: 

Min. 1st Qu. Median Mean 3rd Qu. Max. 

-2.840000 0.000000 0.000000 -0.002114 0.000000 3.686000 

> plotCres.cp, main = "Compound Poisson Increment of a CDGARCH(1,1) model") 


I Increment of a COGARCH(1,1) model 



We are able also to generate the original process using the increments stored into the object 
res.cp using the simulate function. 

> traj.cp<- simulate(res.cp) 

> plot(traj.cp, main = "estimated C0GARCH(1,1) driven by compound poisson process") 

In the next example, we simulate and estimate a COGARCH(l,l) model driven by a 
Variance Gamma process. We set the values for the parameters and dehne the model using 
the following command lines: 

> param.VG <- list(al = 0.038, bl = 0.053, aO = 0.04/0.053, 

+ lambda = 1, alpha = sqrt(2), beta = 0, mu = 0, xOl = 50.33) 

> cog.VG <- setCogarch(p = 1, q = 1, work = FALSE, 

+ measure=list("rngamma(z, lambda, alpha, beta, mu)"). 
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+ measure.type = "code", Cogarch.var = "y", V.var = "v", 

+ Latent.var = "x", XinExpr = TRUE) 

We obtain a trajectory for C0GARCH(1,1) with Variance Gamma noise. 

> set.seed(numSeed) 

> samp.VG <- setSampling(Terminal = Term, n = num) 

> sim.VG <- simulate(cog.VG, true.parameter = param.VG, 

+ sampling = samp.VG, method = "euler") 

> plot(sim.VG, main = "simulated C0GARCH(1,1) model driven by a Variance Gamma process" 



and then we estimate the model parameters: 

> res.VG <- gmm(sim.VG, start = param.VG, Est.Incr = "IncrPar") 

> summary(res.VG) 


Two Stages GMM estimation 
Call: 

gmm(yuima = sim.VG, start = param.VG, Est.Incr = "IncrPar") 


Coefficients: 


Estimate 
bl 0.051449188 
al 0.028791052 
aO 1.248576654 
lambda 1.049274382 
alpha 1.466220182 
beta 0.051526860 
mu 0.003357025 


Std. Error 
0.04168365 
0.01810412 
NA 

0.09432438 
0.08769087 
0.03929050 
0.02935151 
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Log.objFun L2: -3.755496 


Number of increments: 24000 
Average of increments: 0.003635 
Standard Dev. of increments: 0.258127 


-2 log L of increments: 4291.499302 

Summary statistics for increments: 

Min. 1st Qu. Median Mean 3rd Qu. Max. 

-5.548000 -0.001729 0.000000 0.003635 0.002092 4.005000 

> plotCres.VG, main = "Variance Gamma Increment of a C0GARCH(2,1) model") 


Variance Gamma Ineremenl of a C0GARCH(2,1) model 



Even in this case we can obtain the COGARCH(l,l) trajectory using the estimated incre¬ 
ments as follows: 

> traj.VG <- simulate(res.VG) 

> plotCtraj.VG, main="estimated C0GARCH(1,1) model driven by Variance Gamma process" 
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6.2 COGARCH(p,q) model driven by a Compound Poisson 
process 

We conclude this Section by illustrating an example using COGARCH(p,q) process. In this 
way, we show the ability of the yuima package in managing in a complete way these models. 

For this reason in the following we consider a COGARCH(2,l) driven by Gompound Poisson 
Processes where the jump size is normally distributed. 

We define the GOGARGH(2,l) model in the yuima using the command lines: 

> param.cp2 <- listCaO = 0.5, al = 0.1, bl =1.5, b2 = 0.5, 

+ lambda = 1, eta = 0, sig2 = 1, xOl = 2.5, x02 = 0) 

> mod.cp2 <- setCogarchCp = 1, q = 2, work = FALSE, 

+ measure = list(intensity = "lambda",df = list("dnorm(z,eta,sig2)")), 

+ measure.type = "CP", Cogarch.var = "y", V.var = "v", 

+ Latent.var = "x", XinExpr = TRUE) 

We simulate a trajectory. 

> samp.cp2 <- setSampling(Terminal = Term, n = num) 

> set.seed(numSeed) 

> sim.cp2 <- simulate(mod.cp2, sampling = samp.cp2, 

+ true.parameter = param.cp2, method="euler") 

> plot(sim.cp2, main = "simulated C0GARCH(2,1) model driven by a Compound Poisson process") 
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We estimate the model parameters and recover the underlying Levy noise increments: 

> res.cp2 <- gmm(yuima = sim.cp2, start = param.cp2, Est.Incr = "IncrPar") 

> summary(res.cp2) 

Two Stages GMM estimation 
Call: 

gmm(yuima = sim.cp2, start = param.cp2, Est.Incr = "IncrPar") 

Coefficients: 



Estimate 

Std. Error 

b2 

0.0569630413 

0.20054247 

bl 

0.9520642366 

3.54500502 

al 

0.0281299955 

0.09775311 

aO 

0.2956658497 

NA 

lambda 

1.0423762156 

NA 


22 




eta 0.0002425553 
sig2 0.8154399532 


NA 

NA 


Log.objFun L2: -3.323979 

Number of increments: 24000 
Average of increments: -0.001929 
Standard Dev. of increments: 0.258830 

-2 log L of increments: 2861.417140 

Summary statistics for increments: 

Min. 1st Qu. Median Mean 3rd Qu. Max. 

-3.054000 0.000000 0.000000 -0.001929 0.000000 3.483000 

> plot(res.cp2, main = "Compound Poisson Increment of a C0GARCH(2,1) model") 
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The path of the COGARCH(2,l) driven by the estimated increments are reported below: 

> traj.cp2 <- simulate(res.cp2) 

> plot(traj.cp2, main = "estimated C0GARCH(2,1) model driven by a Compound Poisson process") 
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